%	Generate Debt policy functions and asset pricing measure for the Example of Section IX of LANCIA, RUSSO, and WORRALL (2024)   %

clear all; close all;
tic
warning('off')

cd ..\
cd Stored_file

load Canonical_policy

cd ..\
cd Matlab_functions

% autarky levels
eo1=e1-ey1;
eo2=e2-ey2;

v1aut = (utility(ey1,gamma)+betta*(pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma))); % intertemporal utility in autarky in state 1
v2aut = (utility(ey2,gamma)+betta*(pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma))); % intertemporal utility in autarky in state 2
waut = pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma);

% restrict on the relevant state space
wgrid_new = [w_thr:0.001:max(wgrid)];
cy1=interp1(wgrid,cy1,wgrid_new,'spline');
cy2=interp1(wgrid,cy2,wgrid_new,'spline');
w1=interp1(wgrid,w1,wgrid_new,'spline');
w2=interp1(wgrid,w2,wgrid_new,'spline');
wgrid = wgrid_new;

figure('name','policies')
subplot(2,2,1)
plot(wgrid,w1,'r'); hold on;
plot(wgrid,w2,'b'); hold on;
plot(wgrid,wgrid,'black'); hold on;
title('Promised utilities','fontsize',08)
subplot(2,2,2)
plot(wgrid,cy1,'r'); hold on;
plot(wgrid,cy2,'b'); hold on;
title('Individual Consumption','fontsize',08)

% next period consumption
cy11=interp1(wgrid,cy1,w1,'spline');
cy12=interp1(wgrid,cy2,w1,'spline');
cy21=interp1(wgrid,cy1,w2,'spline');
cy22=interp1(wgrid,cy2,w2,'spline');

figure('name','future consumption')
plot(wgrid,cy11,'r'); hold on;
plot(wgrid,cy12,'r--'); hold on;
plot(wgrid,cy21,'b'); hold on;
plot(wgrid,cy22,'b--'); hold on;

% stochastic discount factor
m11=betta.*(cy1)./(1-cy11);
m12=betta.*(cy1)./(1-cy12);
m21=betta.*(cy2)./(1-cy21);
m22=betta.*(cy2)./(1-cy22);

figure('name','SDF')
plot(wgrid,m11,'r'); hold on;
plot(wgrid,m12,'r--'); hold on;
plot(wgrid,m21,'b'); hold on;
plot(wgrid,m22,'b--'); hold on;
title('SDF','fontsize',08)

%%%%%%%%%%%%%%%%%%%%%%%%% FIRST BEST ALLOCATION %%%%%%%%%%%%%%%%%%%%%%%%%%%
% debt as share of young endowment
b1fb=(ey1-delta/(delta+betta))./ey1;
b2fb=(ey2-delta/(delta+betta))./ey2;

cy1fb=delta/(delta+betta);
cy2fb=delta/(delta+betta);

phiy1fb = (cy1fb./ey1);
phiy2fb = (cy2fb./ey2);

phio1fb = ((1-cy1fb)./ey1);
phio2fb = ((1-cy2fb)./ey2);

% bonds' state price 
q11fb = pi1.*betta.*(phiy1fb./phio1fb);
q21fb = pi1.*betta.*(phiy2fb./phio1fb);
q12fb = pi2.*betta.*(phiy1fb./phio2fb);
q22fb = pi2.*betta.*(phiy2fb./phio2fb);

% first best primary surplus/deficit
tau1fb = b1fb - (q11fb.*b1fb+q12fb.*b2fb);
tau2fb = b2fb - (q21fb.*b1fb+q22fb.*b2fb);

% one period-bond price modified=delta
p1fb = q11fb*ey1/ey1+q12fb*ey1/ey2;
p2fb = q21fb*ey2/ey1+q22fb*ey2/ey2;

% one period-bond yield modified
y1fb = -log(p1fb);
y2fb = -log(p2fb);

%%%%%%%%%%%%%%%%%%%%%%%% SECOND BEST ALLOCATION %%%%%%%%%%%%%%%%%%%%%%%%%%%
% critical level of debt
bcrit=1-exp(-betta*(w_thr-waut));

% state-contingent outstanding bond
b1 = (ey1-cy1)./ey1;
b2 = (ey2-cy2)./ey2;

% new issued state-contingent bond
b11 = (ey1-cy11)./ey1;
b12 = (ey2-cy12)./ey2;
b21 = (ey1-cy21)./ey1;
b22 = (ey2-cy22)./ey2;

figure('name','debt')
plot(b1,b11,'b',b1,b12,'r',b1,b1,'k:')

% inverse of marginal utility
phiy1 = (cy1./ey1);
phiy2 = (cy2./ey2);

phio11 = ((1-cy11)./ey1);
phio12 = ((1-cy12)./ey2);
phio21 = ((1-cy21)./ey1);
phio22 = ((1-cy22)./ey2);

phiy11 = ((cy11)./ey1);
phiy12 = ((cy12)./ey2);
phiy21 = ((cy21)./ey1);
phiy22 = ((cy22)./ey2);

% SDF
m11 = betta.*(phiy1./phio11);
m12 = betta.*(phiy1./phio12);
m21 = betta.*(phiy2./phio21);
m22 = betta.*(phiy2./phio22);

% SDF_A: between adjacent generations
m11A = delta.*(phiy1./phiy11);
m12A = delta.*(phiy1./phiy12);
m21A = delta.*(phiy2./phiy21);
m22A = delta.*(phiy2./phiy22);

% SDF_B: between young and old
m11B = (betta/delta).*(phiy11./phio11);
m12B = (betta/delta).*(phiy12./phio12);
m21B = (betta/delta).*(phiy21./phio21);
m22B = (betta/delta).*(phiy22./phio22);

figure('name','Decomposition SDF')
plot(b1,m11,'r',b1,m12,'b')
hold on
plot(b1fb,delta,'r*',b2fb,delta,'b*')
hold on
plot(b1,m11A,'r--',b1,m12A,'b--')
hold on
plot(b1,m11B,'r:',b1,m12B,'b:')
hold on
plot(b1,delta.*ones(1,length(b1)),'k:')

% state price
q11 = pi1.*m11;
q12 = pi2.*m12;
q21 = pi1.*m21;
q22 = pi2.*m22;

% one period-bond price
p = q11+q12;

% bond revenue
BR1=q11.*b11;
BR2=q12.*b12;
TBR=BR1+BR2;
BRmin=interp1(b1,TBR,bcrit,'spline');

% primary surplus
tau = b1 - TBR;

figure('name','BR')
plot(b1,tau,'k',b2,b2,':')
plot(b1, q11.*b11+q12.*b12,'b',b1,b1,'k')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%FIGURES%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure('name','Debt Policy')
% Debt policy
plot(b1,b11,'r','LineWidth',1)
hold on
plot(b1,b12,'b','LineWidth',1)
hold on
plot(b1,b1,'k:')
hold on
plot(bcrit.*ones(1,length(b1)),b1,'k:')
hold on
plot(b2fb.*ones(1,length(b1)),b1,'k:')
hold on
plot(b1fb.*ones(1,length(b1)),b1,'k:')
hold on
plot(min(b2).*ones(1,length(b1)),b1,'k:')
hold on
plot(b1,min(b11).*ones(1,length(b1)),'k:')
hold on
plot(b1fb,b1fb,'r o',b2fb,b2fb,'b o')
hold off
xlim([b1fb max(b11)])
ylim([b1fb max(b11)])
title("Optimal Public Debt Policy")
legend('r=s(1)','r=s(2)')
xlabel('d')

figure('name','Bond Revenue')
% Bond revenue
plot(b1,TBR,'k','LineWidth',1)
hold on
plot(b1,b1,'color', [200 200 200]/255,'LineWidth',1)
hold on
plot(bcrit.*ones(1,length([tau1fb:0.01:max(b1)])),[tau1fb:0.01:max(b1)],'k:')
hold on
plot((bcrit-0.044).*ones(1,length([tau1fb:0.01:max(b1)])),[tau1fb:0.01:max(b1)],'k:')
hold on
plot(b1,(BRmin).*ones(1,length(b1)),'k:')
hold off
xlabel("d")
xlim([b1fb max(b12)])
ylim([min(b1) max(b12)])
title("Bond Revenue")

figure('name','FRF')
% Fiscal reaction function
plot(b1,tau,'k','LineWidth',1)
hold on
plot(b1,0.*ones(1,length(b1)),'k:','LineWidth',1)
hold on
plot(b1fb,tau1fb,'r *',b2fb,tau2fb,'b *')
hold on
plot(bcrit.*ones(1,length([tau1fb:0.01:tau2fb])),[tau1fb:0.01:tau2fb],'k:')
hold on
xlabel("d")
ylabel("\tau(d)")
xlim([b1fb max(b12)])
ylim([tau1fb tau2fb])
title("Fiscal Reaction Function")

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%% MODIFIED STATE SPACE: DEBT AS STATE VARIABLE %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nb = 1000;
wmax=max(wgrid);
bmin = interp1(wgrid,b11,w_thr,'spline'); 
bmax = interp1(wgrid,b22,wmax,'spline');
b1fb = (ey1-delta/(delta+betta))./ey1;
b2fb = (ey2-delta/(delta+betta))./ey2;

% preparation of space
bcrit=1-exp(-betta*(w_thr-waut));
wcrit=fsolve(@(x)interp1(wgrid,b1,x,'spline')-bcrit,-0.7);
wmin = w_thr;
wmax = fsolve(@(x)interp1(wgrid,b12,x,'spline')-b2fb,-0.5); 
wgridmod = [wmin:(wmax-wmin)./(nb-1):wmax];

for i=1:nb
debt(i)=interp1(wgrid,b1,wgridmod(i),'spline');
debt1(i)=interp1(wgrid,b11,wgridmod(i),'spline');
debt2(i)=interp1(wgrid,b12,wgridmod(i),'spline'); 
end
 
% state prices
for i=1:nb
stateq1(i)=betta*pi1*(1-debt(i))./((eo1/ey1)+debt1(i));
stateq2(i)=betta*pi2*(1-debt(i))./((eo2/ey2)+debt2(i)); 
end

figure('name','checks debt')
plot(b1,b11,'r--',b1,b12,'b--');hold on
line([bcrit bcrit],[b1fb b2fb],'linestyle', ':'); hold on
line([b2fb b2fb],[b1fb b2fb],'linestyle', ':'); hold on
plot(debt,debt1,'r',debt,debt2,'b',debt,debt,':')

% risk-adjusted SDF as function of debt
m11b=betta*ey1.*(1-debt)./(1-ey1*(1-debt1));
m12b=betta*ey1.*(1-debt)./(1-ey2*(1-debt2));
m21b=betta*ey2.*(1-debt)./(1-ey1*(1-debt1));
m22b=betta*ey2.*(1-debt)./(1-ey2*(1-debt2));

% variance of SDF
Em1b=pi1.*m11b+pi2.*m12b;
Em2b=pi1.*m21b+pi2.*m22b;

Em1b_2=pi1.*m11b.^2+pi2.*m12b.^2;
Em2b_2=pi1.*m21b.^2+pi2.*m22b.^2;

Varm1b=Em1b_2-Em1b.^2;
Varm2b=Em2b_2-Em2b.^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RISK PREMIUMS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% RISK PREMIUM ON AGGREGATE RISK %
sigma=0.23; % variance of aggregate shock
%mu_gamma=50/49;
mu_gamma=1.04;
pgH=0.5;
pgL=1-pgH;

gH=mu_gamma+sigma;                   % high growth 
gL=mu_gamma-sigma.*(pgH/pgL);        % low growth

Eg=pgH.*gH+pgL.*gL;                  % E(gamma)
Vg=(sigma.^2.*(pgH/(1-pgH)));        % V(gamma)

Eg1=pgH.*gH.^(-1)+pgL.*gL.^(-1);     % E(1/gamma)
Eg2=pgH.*gH.^(-2)+pgL.*gL.^(-2);     % E(1/gamma)
Egg1=Eg.*Eg1;                        % E(gamma)*E(1/gamma)
Vg1=(sigma.^2.*(pgH/(1-pgH))).*(1./((mu_gamma+sigma).*(mu_gamma-sigma.*(pgH)/pgL))).^(2); % V(1/gamma)
harmonic=1./Eg1;
MRP_fb=-(1-Eg.*Eg1);                 % Aggregate MRP Full risk sharing

% (RISK ADJUSTED) MULTIPLICATIVE RISK PREMIUM ON DEBT RETURN
co1d=1-ey1.*(1-debt1);
co2d=1-ey2.*(1-debt2);

ERd1=(1./(betta.*ey1.*(1-debt))).*((pi1.*ey1.*debt1+pi2.*ey2.*debt2)./(pi1.*ey1.*debt1.*(1./co1d)+pi2.*ey2.*debt2.*(1./co2d)));
ERd2=(1./(betta.*ey2.*(1-debt))).*((pi1.*ey1.*debt1+pi2.*ey2.*debt2)./(pi1.*ey1.*debt1.*(1./co1d)+pi2.*ey2.*debt2.*(1./co2d)));

R1f=(1./(betta.*ey1.*(1-debt))).*(1./(pi1.*(1./co1d)+pi2.*(1./co2d)));
MRP=(ERd1-R1f)./R1f;

alpha=ERd1./R1f;
MRP_fe=-(1-Eg.*Eg1);
MRPAgg=MRP+alpha.*MRP_fe;

figure('name','MRP')
subplot(2,2,1)
plot(debt,MRPAgg,'k',debt,MRP,'k--')
title('MRP')
xlabel('d')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% RECURSIVE PRICE OF BONDS %%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nn=10;
p_0(1,:) = ones(1,length(debt));

for j=2:nn
for i=1:length(debt)
p_0(j,i)=stateq1(i).*interp1(debt,p_0(j-1,:),debt1(i),'spline')+stateq2(i).*interp1(debt,p_0(j-1,:),debt2(i),'spline');
end
end

plot([1:1:nn],p_0(:,1),'k',[1:1:nn],p_0(:,length(debt)),'k*',[1:1:nn],0.*ones(1,nn),'k:')

% yield 
for j=2:nn
y0(j,:)=-(1/(j-1)).*log(p_0(j,:));
end                                       

y0=y0([2:1:nn],:);

% one period excess return of bond
for j=2:nn
for i=1:length(debt)
psi1(j,i)=log(interp1(debt,p_0(j-1,:),debt1(i),'spline')/p_0(j,i))+log(delta);
psi2(j,i)=log(interp1(debt,p_0(j-1,:),debt2(i),'spline')/p_0(j,i))+log(delta);

r1(j,i)=(interp1(debt,p_0(j-1,:),debt1(i),'spline')/p_0(j,i));
r2(j,i)=(interp1(debt,p_0(j-1,:),debt2(i),'spline')/p_0(j,i));
end
end

figure('name','yields')
% yield curve
subplot(2,2,1)
plot([1:1:nn-1],-log(delta)*ones(1,nn-1),'linestyle', ':', 'color', 'k','LineWidth',0.1)
hold on
plot([1:1:nn-1],y0([1:1:nn-1],1),'Color',[0.8 0.8 0.8],'LineWidth',1.5)
hold on
plot([1:1:nn-1],y0([1:1:nn-1],60),'k--','LineWidth',1)
hold on
plot([1:1:nn-1],y0([1:1:nn-1],length(debt)),'k','LineWidth',1)
hold off
legend('-log(\delta)','d=d_{min}','d_{min}<d<b^{\ast}(2)','d=b^{\ast}(2)')
ylabel('i^k(d)')
xlabel('k')

% one period long bond excess returns
subplot(2,2,2)
plot([1:1:nn-1],0*ones(1,nn-1),'linestyle', ':', 'color', 'k','LineWidth',0.1)
hold on
plot([1:1:nn-1],psi2([1:1:nn-1],[1]),'Color',[0.8 0.8 0.8],'LineWidth',1)
hold on
plot([1:1:nn-1],psi2([1:1:nn-1],[length(debt)]),'--','Color',[0.8 0.8 0.8],'LineWidth',1.5)
hold on
plot([1:1:nn-1],psi1([1:1:nn-1],[1]),'k','LineWidth',1)
hold on
plot([1:1:nn-1],psi1([1:1:nn-1],[length(debt)]),'k--','LineWidth',1.5)
hold on
plot([1:1:nn-1],0.7142.*ones(1,nn-1),'r--','LineWidth',0.5)
hold on
plot([1:1:nn-1],-0.7142.*ones(1,nn-1),'r--','LineWidth',0.5)
hold off
legend('0','\psi^k(b^{\ast}(1),2)','\psi^k(b^{\ast}(2),2)','\psi^k(b^{\ast}(1),1)','\psi^k(b^{\ast}(2),1)')
ylabel('\psi^k(d,s^{\prime})')
xlabel('k')

% inverse of return of long bond
subplot(2,2,3)
plot(debt,r1(nn,:).^(-1),'r')
hold on
plot(debt,r2(nn,:).^(-1),'b')
hold off

% SDF
subplot(2,2,4)
plot(debt,psi1(nn,:),'r')
hold on
plot(debt,psi2(nn,:),'b')
hold on
plot(debt,-log(stateq1./pi1)+log(delta),'r--')
hold on
plot(debt,-log(stateq2./pi2)+log(delta),'b--')
hold on

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LONG RUN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% TRANSITION MATRIX
edge= debt(1:nb);
edge(nb+1)=debt(nb);
debtk = discretize(debt,edge);
debt1k = discretize(debt1,edge);
debt2k = discretize(debt2,edge);

debt1k(1)=1;

% correct
for i=2:nb
if isnan(debt1k(i))
   debt1k(i)=debt1k(i-1); 
end
if isnan(debt2k(i))
   debt2k(i)=debt2k(i-1);  
end
end

% matrix of transition probabilities and state prices
for i=1:nb
   for j=1:nb
       if debtk(j)==debt1k(i)
            T11(i,j)=pi1;
            Q11(i,j)=stateq1(i);
       else
            T11(i,j)=0;
            Q11(i,j)=0;
       end
   
      if debtk(j)==debt2k(i)
           T12(i,j)=1-pi1;
           Q12(i,j)=stateq2(i);
      else
           T12(i,j)=0;
           Q12(i,j)=0;
      end
      
      if debtk(j)==debt1k(i)
           T21(i,j)=pi1;
           Q21(i,j)=stateq1(i);
      else
           T21(i,j)=0;
           Q21(i,j)=0;
      end
      
      if debtk(j)==debt2k(i)
          T22(i,j)=1-pi1;
          Q22(i,j)=stateq2(i);
      else
          T22(i,j)=0;
          Q22(i,j)=0;
      end
      
   end
end

T=[T11 T12; T21 T22]; 
Q=[Q11 Q12; Q21 Q22];
T=T';

% stationary distribution
probb = (1/(2*nb))*ones(2*nb,1);
   test=1;
   while test > 12^(-12);
      probbx = T*probb;
      test = max(abs(probbx-probb));
      probb = probbx;
   end;

probb=probb/sum(probb); 

for i=1:nb
  probb1(i)=probb(i);
end

sum(probb1)

for i=nb+1:2*nb
   probb2(i-nb)=probb(i); 
end
sum(probb2)

probb = probb1+probb2;

figure('name','LR distribution')
plot(debt,probb,'LineWidth',1.5);
xlim([b1fb max(b12)])
ylim([0 0.25])
set(gca,'XTick',[], 'YTick', [])
text(max(b22)-0.01,0-0.015,'$$d_{\max}$$','Interpreter','Latex','FontSize',12)
text(min(b11)-0.03,0-0.015,'$$d_{\min}=b^{\ast}(1)$$','Interpreter','Latex','FontSize',12)
text(b2fb-0.01,0-0.015,'$$b^{\ast}(2)$$','Interpreter','Latex','FontSize',12)
text(bcrit-0.01,0-0.015,'$$d^c$$','Interpreter','Latex','FontSize',12)
text(min(b11)-0.04,0.125,'$$0.125$$','Interpreter','Latex','FontSize',12)
text(min(b11)-0.04,0.25,'$$0.25$$','Interpreter','Latex','FontSize',12)
text(min(b11)-0.045,0.0625,'$$0.0625$$','Interpreter','Latex','FontSize',12)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%% ROSS MEASURE OF MAX VARIABILITY OF YIELDS %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x0 = ones(length(Q),1);

distmia=100;
iter=0;
tollv=10^(-10);

while distmia>=tollv
iter = iter+1;

x1 = Q*x0;
yy0 = x0./max(x0);
yy1 = x1./max(x1);

distmia = max(abs(((yy1)-(yy0))./((yy0))));

x0 = x1;
end

max_eig = yy1([1:1:length(debt)]);
max_eig1=interp1(debt,max_eig,debt1,'spline');
max_eig2=interp1(debt,max_eig,debt2,'spline');

% Perron root
rho  =(Q*yy1)'*yy1/((yy1)'*yy1);

Ross = -log(min(max_eig)/max(max_eig));

% closed form solution of Ross measure reported in Proposition 8
Gamma = -((1-ey1)/(ey1))+exp((1/(betta*(1-pi2)))*log((ey2*(((betta+delta))/delta))*((1-ey2)*(((betta+delta))/betta))^(betta*pi2)...
    .*(((1-ey1)/(ey1)))^(betta*(1-pi2))));

Ross_geom = log(((ey2)/(ey1))*((Gamma+((1-ey1)/(ey1)))/((betta/delta)*(1-Gamma))));

cd ..\
cd Stored_file

save ('Debt_policy_2state.mat','debt','probb1', 'probb2','MRP_fe','alpha','MRP','b1fb','b2fb')

